Lighting the Universe with filaments 
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The first stars in the Universe form when chemically pristine gas heats as it 
falls into dark matter potential wells, cools radiatively due to the formation 
of molecular hydrogen, and becomes self-gravitating. We demonstrate with 
super-computer simulations that their properties depend critically on the cur- 
rently unknown nature of the dark matter. If the dark matter particles have in- 
trinsic velocities that wipe-out small-scale structure, then the first stars form in 
filaments with lengths of order the free-streaming scale, which can be ~ 10 2() m 
(~ 3 kpc, baryonic masses ~ 10 7 solar masses) for realistic 'warm dark mat- 
ter' candidates. Fragmentation of the filaments forms stars with a range of 
masses which may explain the observed peculiar element abundance pattern 
of extremely metal-poor stars, while coalescence of fragments and stars dur- 
ing the filament's ultimate collapse may seed the super massive black holes 
that lurk in the centres of most massive galaxies. 

Most of the matter in the Universe does not interact with light except gravitationally. This 
'dark matter' is usually assumed to be 'cold', meaning that its velocity dispersion is sufficiently 
small for density perturbations imprinted in the early Universe to persist up to very small scales. 
Although this model is able to describe the large-scale distribution of galaxies in impressive 
detail, it may face problems on the scale of galaxies and below, for example it may predict too 
many satellite galaxies ©, as well as too cuspy profiles for the dark matter halos that surround 
galaxies ©. 

However, dark matter has yet to be detected in the laboratory, and there exist many viable 
dark matter candidates from particle physics that are not 'cold' . Warm dark matter (WDM) 
particles have intrinsic thermal velocities and these motions quench the growth of structure be- 
low a 'free-streaming' scale (the distance over which a typical WDM particle travels) which 
depends on the nature of the particle. Because small and dense halos do not form below the 
free-streaming- scale, the dark matter halos that surround galaxies in a WDM model have far 
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less substructure and are less concentrated as compared to their cold dark matter (CDM) coun- 
terparts, which may help alleviate both the satellite and galactic core problem (3). Structures on 
larger scales are similar in WDM and CDM and therefore the distribution of galaxies is not af- 
fected. The first generation of stars in the Universe forms when primordial gas gets compressed 
by falling into small dark matter potential wells (@|5]|^0). Because WDM affects structure 
formation on such small scales it may influence how the first stars form; we have performed 
simulations to analyse this in more detail. 

Large-scale power in the spectrum of density perturbations causes progenitors of present- 
day clusters of galaxies to be amongst the first objects to condense out of the initially almost 
smooth mass distribution. We study the early formation stages of such an object by identifying 
a massive cluster of galaxies in a dark matter simulation of a large cosmological volume at red- 
shift z = 0, and use a multi-scale technique d#]|9]) to re-simulate its formation and evolution with 
the cosmological hydrodynamics code Gadget-2 (1701191). Baryons compressed by falling into the 
developing dark matter potential wells cool radiatively through molecular hydrogen emission 
lines (|771 IPb: we follow the formation of molecular hydrogen and its cooling using a sophisti- 
cated chemistry network ©. The parent simulation, from which we picked the re-simulated 
region, followed the growth of structure in a CDM dominated Universe with a cosmological 
constant A (ACDM); we performed the re- simulations assuming both CDM and WDM as dark 
matter. In this WDM model, small scale power in the density fluctuations is exponentially sup- 
pressed below ~ 100 co-moving kpc, mimicking free- streaming of gravitino particles with mass 
mwDM = 3 keV ©. Gravitinos are a popular WDM candidate (12) but our results are insensitive 
to the exact choice of WDM particle as long as its free- streaming length is more than a few tens 
of co-moving kilo parsecs. Therefore even if the gravitino were slightly more massive, or if the 
WDM were instead a sterile neutrino (another popular WDM candidate (1731)), our results would 
not change appreciably. Observations of the clustering of neutral gas along sight lines to distant 
quasars (the Lyman-or forest of absorption lines) probe scales 1-40 Mpc when the density per- 
turbations on these scales are still small, and the presence or absence of significant substructure 
in these observations can constrain the masses of WDM particles <\14\\15\i : our choice of WDM 
particle mass (3 keV, the corresponding mass-scale is ~ 3 x 10 8 M o ) is well above this lower 
limit (~ 2 keV; <\14\\15\i ). The initial amplitude of the imposed density perturbations in our sim- 
ulations is normalised to the level seen in the Cosmic Microwave Background radiation (1731) , 
and our simulations start at redshift z ~ 200. (See the supplementary online material for more 
details on these simulations.) 

The growth of structure in these re-simulations leads to a pattern of filaments and sheets 
(Fig.l) familiar from the local large-scale distribution of galaxies. This is because the assumed 
Gaussian spectrum of density perturbations, appropriate for an inflationary model, leads to col- 
lapse along one (sheet) and two directions (filament) before the formation of halos. Although 
the large-scale filamentary pattern is very similar in CDM and WDM, the structure of the fila- 
ments themselves is very different: whereas the CDM filaments fragment into numerous nearly 
spherical high density regions ('halos'; panel l.a), the WDM filaments are mostly devoid of 
such substructure (panel l.b). Panels l.c and l.d depict the WDM filament at an earlier time 
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before any of the other filaments have formed: the central density is very high (hydrogen num- 
ber density n h ~ 10 4 cm" 3 « 10 6 (/i h )) yet no dark matter halo has formed yet. It is well known 
that the Poisson noise in simulation codes that use particles to represent the dark matter leads 
to spurious fragmentation of the filaments that form in such WDM simulations (|771[7^l) . We 
therefore end the analysis of our WDM simulations well before filaments fragment. 

The length of the filament (~ 3 kpc) is of order of the imposed WDM free-streaming scale 
as expected, L ~ 4 kpc at redshift z = 23.34 when the Universe is 140 million years old. Gas 
and dark matter accrete perpendicular onto the filament's axis (Fig. 2). Dark matter particles 
falling into the filament perform damped oscillations as the potential well deepens. At r ~ 50 pc 
(where r is the distance perpendicular to the filament's axis) dark matter particles falling into the 
well encounter particles that fell in from the other side, and such successive instances of 'orbit- 
crossing' give rise to the steps in the density seen in Fig. lb. Baryons do not undergo orbit- 
crossing but the gas gets compressed to a temperature T x 7000 K at r ~ 20 pc. Rapid build-up 
of H 2 induces cooling and the gas starts to dominate the matter density further downstream, so 
that the ratio of gas to dark matter densities PblPmA ~ 15 = 100(pt,)/(p DM ) at r = 2 pc. At 
r < 2 pc where the gas dominates, the ratios of principal axes of the filament are b/a = 0.123 
and c/a = 0.118 hence the filament is very nearly cylindrical. The properties of both the gas 
and the dark matter are very uniform along the whole length of the filament. The cylindrical 
density profile below 10 pc is approximately p oc r -2,8 for 2 < r < 8 pc and p oc r~ 2 for r < 2 pc. 
This contrasts with p oc r~ 23 for the spherically averaged profile of the gas in CDM halos on a 
comparable scale (00). The central H 2 abundance reaches 10" 3 , higher than in the CDM case 
because of the higher temperature reached behind the accretion shock. This higher temperature, 
and the associated higher ionisation fraction, will also enhance the importance of HD cooling 
at later stages (li 91 1201) . 

The non-linear collapse into a thin filament found in these WDM simulations is in sharp 
contrast to what happens in the CDM case. There the first objects to reach high densities 
are discrete, nearly spherical, dark matter halos (gravitationally bound concentrations of dark 
matter) that form at tiny masses and build-up hierarchically through mergers and accretion. 
Some halos have a sufficiently deep potential well to accrete and shock baryons, enabling H 2 
formation and radiative cooling. Runaway collapse of the rapidly accreting self-gravitating gas 
is thought to lead to the formation of a single massive star per cooling halo (000). The absence 
of small scale power in WDM prevents halos from forming before the filament itself becomes 
highly non-linear. Although our re-simulations focus on the progenitor of a massive cluster, 
and hence the forming object collapses unusually early on, the fact that very high densities are 
reached in a filament, as opposed to a spherical halo, is generic to WDM. 

The stability of collapsing filamentary clouds has been investigated in the context of the 
formation of cloud cores by for example (|271l2"2)) . and applied to early Universe filaments (1231) . 
The inability of gas to cool sufficiently fast usually limits the collapse time of the filament by 
the cooling time, td = pi 'p ~ t c = nk Bo \ tz T/A (where T is the temperature and A is the cooling 
rate due to H 2 and HD cooling). Perturbations start growing, possibly leading to fragmentation, 
when the dynamical time t d » t p , where t p is the inverse growth rate of perturbations. If gas 
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cooling is efficient, td ~ t c «: t p , and perturbations do not grow. At densities n = n\ ~ 10 5 cm -3 
the level population of H 2 reaches LTE making A oc n instead of oc n 2 and t d increases (|77T) . 
A sufficiently massive filament may yet survive fragmentation at this stage, and at higher den- 
sity n ~ 10 9 cm" 3 three-body processes promote the formation of H 2 and the cooling time 
decreases again. However when n > rc 2 = 10 12 cm -3 the gas becomes optically thick in the 
H 2 cooling transitions slowing the collapse and the filament is once more in danger of frag- 
menting. Collision induced continuum emission will again decrease the cooling time when 
n ~ 10 14 cm" 3 , until the gas becomes optically thick also to this cooling radiation at densities 
n > h 3 = 10 16 cm -3 (|2?ll2Ql) . The physics of the H 2 molecule therefore sets three densities 
at which the filament may fragment. The typical fragment masses are of order tens of solar 
masses, solar masses, and sub-solar masses, for fragmentation at densities n\, n 2 and n 3 , respec- 
tively (1231) . 

The tidal field around the filament breaks the cylindrical symmetry on scales comparable to 
the filament's length (Fig. 1), and can trigger the gas dynamical instabilities that ultimately lead 
to fragmentation. Unfortunately our current simulations are not able to follow this process in 
detail because these tiny deviations from symmetry are overpowered by numerical noise caused 
by the graininess of the particle distribution. This causes the filament to fragment very rapidly 
as expected, yet the scale of the fragmentation is artificial. In the WDM universe, the small- 
scale perturbations that trigger fragmentation in the simulations are not present and need to be 
generated through transfer of power from larger scales. Since this is a relatively slow process, 
the central density will have reached the higher value n 2 or even higher, implying small fragment 
masses of order of a solar mass or below. Such fragments can coalesce to form more massive 
clumps, as demonstrated in 2D numerical simulations (|22l) . Even if individual cores survive, 
they may still grow in mass through accretion. 

Detailed observations of star-forming clouds in the Milky Way reveal that the low-mass stel- 
lar mass function is very similar to that of the dense, pre-stellar cloud cores within the ambient 
cloud, although it is not yet clear what determines the cut-off at high masses. This close similar- 
ity might indicate that the process of cloud fragmentation plays an important role in determining 
the initial mass function (1251) . If this also applies to star-formation in a WDM filament, then it 
is plausible that fragmentation will lead to a burst of star formation which includes low-mass 
stars with masses ~ 1 M Q or below, but also much more massive stars built through mergers 
and accretion. The more massive stars will end their short lives through supernova explosions 
or collapse to form intermediate mass black holes, and some of the very low mass stars may 
potentially survive until today. Although the details of this scenario are uncertain, it is clear that 
the stellar mass function will be quite different from the CDM case. 

The low-mass Milky Way stars HE 0107-5240 (26) and HE 1327-2326 (^3 have extremely 

low metallicities^ of [Fe/H] 5 and peculiar element abundances, for example [C/Fe]>l. An 

initial mass function of zero metallicity first stars with a range of masses, as we suggested might 
be the case in WDM, could explain these stars either as due to self-pollution or pre-enrichment 

'The notation [X] = log(X/Xo), characterises the mass fraction of element X in terms of the value in the Sun, 

X . 



4 



by intermediate mass (tens of solar masses) supernovae (|2ffl) . Such first generation low and 
intermediate mass stars are not thought to form in CDM (|4]|6][Z]>- Therefore, although the 
present observational evidence is not yet unambiguous, a future detection of zero-metallicity 
low-mass first stars may indicate that the dark matter is warm. 

Free-streaming of the WDM decreases significantly the number density of halos that host 
early star formation, which could delay reionization as compared to CDM (29). However the ad- 
ditional mode of star formation in filaments could partly compensate, making it unclear whether 
reionization is indeed delayed when the dark matter is warm. 

What is the ultimate fate of such ~ 10 7 M Q filaments? Eventually the filament will collapse 
along its long axis, and since its mean density is very high a significant number of collisions 
between cloud cores and stars would appear inevitable. Such collisions could build-up a massive 
object which can seed the formation of the super-massive black holes that power redshift z ~ 6 
quasars, and appear to lurk in the centers of most large galaxies today. 

The different outcome of WDM versus CDM is because the thermal velocities of the WDM 
particles prevent halos from forming before the filaments themselves form stars. This happens 
when the free-streaming length introduced by the dark matter's thermal velocities is more than a 
few tens of co-moving kilo parsecs so that the filaments are sufficiently massive to heat infalling 
gas to more than ~1000 K, enabling efficient cooling by molecular hydrogen. If we simulate a 
WDM Universe with much shorter free- streaming length, for example ~ 20 kpc corresponding 
to a gravitino mass m W DM = 15 keV, star formation proceeds in a similar way to the CDM case, 
in good agreement with earlier work (30). The likely very different initial mass function of the 
first stars and the rapid formation of massive black holes in a WDM scenario with m WDM ~ 
3 keV, as opposed to CDM, implies a very different early thermal and metal enrichment history 
of the Universe, greatly affecting subsequent galaxy formation. It appears therefore that the way 
in which quasar, star and galaxy formation started, depends strongly on the nature of the dark 
matter. 
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Figure 1 

Top panels: Dark matter density structure of the progenitor of a redshift z = massive cluster 
of galaxies at z = 15 when the thermal velocities of the dark matter particles are negligible 
(CDM; panel a) and when the dark matter is 'warm' (WDM; a gravitino with m W DM = 3 keV; 
panel b). Although both models produce a characteristic filamentary pattern in the density, the 
CDM filaments fragment into numerous nearly spherical halos whereas free- streaming of the 
WDM prevents such substructure from forming. Bottom panels: Gas temperature (panel c) and 
density (panel d) in the WDM filament indicated by the red box in panel b, at an earlier redshift 
(z = 23.34) when only this filament had formed. Gas accretes very uniformly onto the filament 
as indicated by the velocity vectors, heats as it gets compressed, but further downstream cools 
due to the formation of H2, making the centre of the filament cold and dense. The filament 
shown here is almost perfectly cylindrical, more generally they are elliptical in shape. 
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Figure 2 

Gas and dark matter profiles of the WDM filament shown in panels c and d of Fig. 1 . Top left: 
accretion of dark matter particles (black) and gas particles (red); note the different horizontal 
and vertical scales (along and perpendicular to the filament, respectively). The dark matter 
undergoes orbit-crossing (multi-valued velocities) at r ~ 50 pc, whereas the gas gets shocked 
much closer in (~ 20 pc) and dominates further downstream. Top right: gas density profile as 
measured perpendicular to the filament's axis. Lines p oc r~ 2 and p oc r" 2 8 are drawn to guide 
the eye. Bottom left: temperature profile: gas heats as it gets compressed, but the rapid build-up 
of H 2 cools the gas further in. Bottom right: Build-up of molecular hydrogen fraction, f Hl , 
during the accretion. 
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We briefly describe the setup of the initial conditions, how we model the effect of warm dark 
matter, the chemistry network, and the numerical method used. 

1 Numerical Methods 

1.1 Cosmological initial conditions 

The theory of inflation posits that quantum mechanical fluctuations in the 'inflaton' field were 
blown-up to macroscopic scales during a short period of rapid expansion called inflation. These 
small perturbations in the density have grown due to gravity into the large-scale structure we see 
today (SI). The initial conditions of a cosmological simulation represent these initially small 
perturbations at the starting redshift of the simulation in a periodic, computational box of given 
co-moving size. The simulation then follows how these small ripples develop into high and 
low-density regions. 

If the statistics of the perturbations are those of a Gaussian field, as is thought to be the case 
in inflation, then the initial fluctuations are fully described by their power- spectrum, P(k, t), 
which characterises the amplitude of the fluctuations as function of scale, A = In Ik, at time 
t. The Warm Dark Matter (WDM) and Cold Dark Matter (CDM) power spectra from (S2) are 
plotted in Fig. SI: free- streaming of a WDM particle with mass m W DM = 3 keV suppresses 
power below a scale ~ 100 kpc as compared to CDM. 

Setting-up initial conditions now boils down to generating a Gaussian field with known 
power-spectrum, which is done using Fourier transform methods. Since our simulation code 
uses particles to represent the matter distribution, we need to generate a particle distribution to 
represent this density field. We start by creating a homogeneous particle distribution, either by 
simply putting equal mass particles on the vertices of a regular grid, or as done here, starting 
from a glass distribution. We then use the Zel'dovich approximation (S3) to displace these 
particles to represent the density field (see also S4). The physical power- spectrum has power 
on all scales, but the particle distribution cannot represent waves larger than the box size, or 
smaller than twice the inter particle spacing d. 

The dark matter structures that host the very first stars are very rare objects, therefore cos- 
mological simulations of such early structure formation must model extremely large volumes 
in order to sample these rare objects. However, the simulations must also achieve very high 
mass resolution to follow accurately how gas collapses to high densities in these early potential 
wells. The contradictory requirements of large simulation volume yet high resolution make cos- 
mological simulation of these structures unfeasible with uniform mass resolution given current 
computational resources. Therefore we use a 'multi-scale re-simulation technique' to overcome 
this dilemma. In a re-simulation, only a sub-volume of the whole simulation volume (contain- 
ing the object in which we are interested) is simulated with very high resolution, while the rest 
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of the computational domain is sampled at much coarser resolution. The advantage is that it en- 
ables us to obtain the extremely high dynamic range in mass and length required for the current 
problem; the challenge is how to generate the appropriate initial conditions (ICs). 

We begin by generating ICs in a large periodic box (479 Mpc, the 'parent simulation') as 
described above, and evolve the system forward in time using Gadget-2 (see S5), following just 
the growth of structure in the dark matter for simplicity and speed. Once this simulation reaches 
redshift z = 0, we analyse it and search for a massive dark matter halo, such as would envelope 
a massive cluster of galaxies (or any other structure of interest). 

Next we create a multi-mass particle distribution to represent the unperturbed homogeneous 
universe, by concentrating lower-mass particles in the region of interest, with successive layers 
of more massive particles away from this region. The original displacements field for the parent 
simulation is recreated and the displacements applied everywhere. In the region of interest, the 
inter particle spacing is smaller, and we can include power on smaller scales than was possible 
with the original coarser particle distribution. This small scale power is generated by creating a 
periodic cubic Fourier mesh which is placed around the region of interest. A displacement field 
is created on this mesh, with the appropriate power spectrum, taking wave numbers between 
the cut-off in power spectrum of the parent simulation and the Nyquist frequency, k = n/d, of 
the region of interest determined by the inter particle spacing, d, there. Typically the longest 
waves used in the small scale power have a maximum wavelength of a tenth of the cubic Fourier 
mesh (if the maximum wavelength were to be much larger this procedure would start becoming 
inaccurate because periodicity is enforced by the method over the mesh). The second set of 
displacements which represent the additional small scale power, is then applied to the particles 
inside the high resolution region only. The scheme used to interpolate the displacements to the 
positions of the particles takes account of the mass of the particles to avoid aliasing effects for 
the most massive particles in the high resolution region. 

Having created these initial conditions and run them forward in time, a new region of interest 
can be defined for this re-simulation at the same or different redshift, and a re-simulation of the 
re-simulation can be made following the same procedure where the parent simulation is now 
identified as the first re-simulation. The simulations described in this paper apply this procedure 
four times in order to generate ICs for the progenitor of a massive cluster of galaxies at redshift 
z = 0, and its surroundings; see (S6) for further details. The main difference here is that we 
used a WDM power spectrum for generating the initial density field. Note, for simplicity we 
did not add the thermal velocities of the WDM particles, hence the free-streaming of the WDM 
is mimicked just by their suppression of small-scale power. 

The highest resolution achieved in our final simulation is M dm = 272.6 M Q and M gas = 
41.9 M Q for dark matter and gas particles respectively, and the high resolution region has a 
Lagrangian radius of about 600 comoving kpc. 
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1.2 The simulation code Gadget-2 



Gadget-2 (see S5) is a Lagrangian code that uses two sets of particles to represent dark mat- 
ter and gas. Each gas and dark matter particle represents a swarm of baryonic particles and 
gravitinos, respectively. Because the simulation particles are so much more massive than the 
physical particles they represent, we soften the gravitational forces between simulation parti- 
cles to minimise numerical artifacts. Each dark matter particle in the simulation has a position, 
velocity, (constant) mass, and gravitational softening length. Gas particles have in addition 
a thermal energy and gas smoothing length, in addition to properties to follow the molecu- 
lar hydrogen formation (see next section). Gadget-2 uses the smoothed particle hydrodynamics 
(SPH) scheme to evaluate pressure gradients between particles. Once gravitational and pressure 
accelerations are computed, the state of system is marched forward in time. Gadget-2 uses a 
sophisticated tree structure to perform these force calculations efficiently on a parallel computer 
with message-passing interface (MPI). Full details on Gadget-2 are given in S5. 



2 Primordial gas chemistry 

In this section, we briefly summarise the basic gas processes that are important for primordial 
gas unpolluted by metals; full details are given in (S6-10). 

At high redshifts, in the absence of metals, and when the temperature of gas is lower than 
10 4 K (above which atomic hydrogen line cooling is dominant), the main coolant enabling star 
formation in the early dark matter potential wells, is molecular hydrogen, H 2 . Its formation has 
three main channels: 

1 . the channel 

H + + H -» H+ + y 

H+ + H -» H 2 + H + 

2. the H channel 

H + e~ -> H + y 

FT + H -» H 2 + e 

3. the three-body reaction 

H + H + H H 2 + H, 

H + H + H 2 -» H 2 + H 2 . 

The first channel dominates at high redshift, z > 200, and the three-body reactions are important 
only at high densities (n H ~ 10 8 cm 3 ); therefore the H channel is most relevant for the initial 
stages in the formation of the filament discussed in this paper. Our simulation code follows 
reactions for 9 species (e~, H, H + , He, He + , He ++ , H 2 , H^, FT), and includes all three H 2 
formation channels. 
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3 Three-dimensional representation of the first filament 



Figure S2 contains projections of the filament shown in Fig. 1 (panels C and D) along all three 
coordinate directions to illustrate its three-dimensional cylindrical form. The structure of the 
filament on larger scales is illustrated in figure S3, where we have transformed the Cartesian 
axes such that the filament lies along the x-axis. The figure illustrates how straight the filament 
is on scales of 100s of parsecs, whereas the tidal field begins to deform the filament on scales 
approaching the free-streaming length. 
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Figure SI: Power- spectrum for CDM (drawn line) and WDM (dashed line). Free- streaming of 
the warm dark matter particle (a gravitino with mass m W DM = 3 keV) exponentially suppresses 
power below a scale ~ 100 kpc. 
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Figure S2: 

Gas temperature (left hand panels) and density (right hand panels) of the filament shown in 
Fig. 1 (panels C and D) along three orthogonal projections. Each panel is 800 (physical) parsec 
on a side. 
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Figure S3: 

Gas temperature (left hand panels) and density (right hand panels) of the filament shown in 
Fig. 1 (panels C and D), after rotating the filament so that its long axis lies along the x-axis. 
Panels from top to bottom zoom out from 4 kr>c to 400 dc. The accretion velocities of the sas 
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